drawing
Matthew Abbott / The New York Times


Introdução

Os incêndios florestais na Austrália são uma ocorrência regular que teve um papel significativo na formação da natureza do continente ao longo de milhões de anos. No entanto, os incêndios podem causar danos materiais significativos e levar à perda de vidas humanas e animais. Os incêndios florestais mataram cerca de 800 pessoas e bilhões de animais desde 1851. Estima-se que entre 2019 e 2020 os incêndios mataram pelo menos 33 pessoas e mais de 3 bilhões de animais.

Geralmente, os incêndios mais devastadores são precedidos por altas temperaturas, baixa umidade relativa e fortes ventos, que criam as condições ideais para a rápida propagação do fogo. A chuva é a maior aliada ao combate ao fogo e na mitigação do surgimento de novos de focos de incêndios. Outro problema que agrava a situação na Austrália, é a forte dependência do país de civis para conter seus incêndios, cerca de 90% dos bombeiros combatendo fogos florestais são voluntários.

Portanto, considerando os recursos limitados para combate aos incêndios florestais, ainda mais com o aumento das mudanças climáticas que estão tornando esses eventos mais frequentes e devastadores, saber se irá chover em uma região específica dadas algumas medições é importante para direcionar os esforços onde eles são mais necessários.

Visualizando os Dados

dataset = tibble(read.csv('weatherAUS.csv'))
print(dataset)
## # A tibble: 145,460 x 23
##    Date       Location MinTemp MaxTemp Rainfall Evaporation Sunshine WindGustDir
##    <chr>      <chr>      <dbl>   <dbl>    <dbl>       <dbl>    <dbl> <chr>      
##  1 2008-12-01 Albury      13.4    22.9      0.6          NA       NA W          
##  2 2008-12-02 Albury       7.4    25.1      0            NA       NA WNW        
##  3 2008-12-03 Albury      12.9    25.7      0            NA       NA WSW        
##  4 2008-12-04 Albury       9.2    28        0            NA       NA NE         
##  5 2008-12-05 Albury      17.5    32.3      1            NA       NA W          
##  6 2008-12-06 Albury      14.6    29.7      0.2          NA       NA WNW        
##  7 2008-12-07 Albury      14.3    25        0            NA       NA W          
##  8 2008-12-08 Albury       7.7    26.7      0            NA       NA W          
##  9 2008-12-09 Albury       9.7    31.9      0            NA       NA NNW        
## 10 2008-12-10 Albury      13.1    30.1      1.4          NA       NA W          
## # ... with 145,450 more rows, and 15 more variables: WindGustSpeed <int>,
## #   WindDir9am <chr>, WindDir3pm <chr>, WindSpeed9am <int>, WindSpeed3pm <int>,
## #   Humidity9am <int>, Humidity3pm <int>, Pressure9am <dbl>, Pressure3pm <dbl>,
## #   Cloud9am <int>, Cloud3pm <int>, Temp9am <dbl>, Temp3pm <dbl>,
## #   RainToday <chr>, RainTomorrow <chr>

Mapa com as Temperaturas Máximas de Cada Local

Usando a API do Google Maps para salvar as coordenadas geográficas de cada local de medição, como existem vários lugares com o mesmo nome em diferentes países anglófonos, é necessário concatenar o nome do país no final do local para garantir que as coordenadas estarão coerentes.

MTemp <- dataset[!is.na(dataset$MaxTemp), ] %>%
  select(Location, MaxTemp)

TempLocation <- aggregate(MTemp$MaxTemp, by=list(Category=MTemp$Location), FUN=mean)
TempLocation$Category <- paste(TempLocation$Category, ", Australia")# Adding Australia to the end of location name since there are many cities with the same names in other countries
cities_df <- as.data.frame(TempLocation)
locations_df <- mutate_geocode(cities_df, Category)
locations <- as_tibble(locations_df)

locations_sf <- st_as_sf(locations, coords = c("lon", "lat"), crs = 4326)
aus_coord <- geocode("Australia")
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Australia&key=xxx
map <- get_googlemap(center = c(aus_coord$lon, aus_coord$lat), zoom = 4)
## Source : https://maps.googleapis.com/maps/api/staticmap?center=-25.274398,133.775136&zoom=4&size=640x640&scale=2&maptype=terrain&key=xxx
ggmap(map) +
  geom_point(data = locations,
             aes(x = lon, y = lat, size = x),
             color = "red", alpha = 0.3)

Mapa de Calor de Pluvialidade

RF <- dataset[!is.na(dataset$Rainfall), ] %>%
  select(Location, Rainfall)

RainLocation <- aggregate(RF$Rainfall, by=list(Category=RF$Location), FUN=sum)
RainLocation$Category <- paste(RainLocation$Category, ", Australia")# Adding Australia to the end of location name since there are many cities with the same names in other countries
Rcities_df <- as.data.frame(RainLocation)
Rlocations_df <- mutate_geocode(Rcities_df, Category)
Rlocations <- as_tibble(Rlocations_df)

Rlocations_sf <- st_as_sf(Rlocations, coords = c("lon", "lat"), crs = 4326)
print(Rlocations_sf)
## Simple feature collection with 49 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 115.1003 ymin: -42.88261 xmax: 167.9547 ymax: -12.46373
## Geodetic CRS:  WGS 84
## # A tibble: 49 x 3
##    Category                       x             geometry
##  * <chr>                      <dbl>          <POINT [°]>
##  1 Adelaide , Australia       4842.  (138.6007 -34.9285)
##  2 Albany , Australia         6828. (117.8837 -35.02693)
##  3 Albury , Australia         5763. (146.9095 -36.07512)
##  4 AliceSprings , Australia   2677. (133.8807 -23.69804)
##  5 BadgerysCreek , Australia  6421.  (150.7634 -33.8751)
##  6 Ballarat , Australia       5269. (143.8503 -37.56216)
##  7 Bendigo , Australia        4913. (144.2794 -36.75702)
##  8 Brisbane , Australia       9941   (153.026 -27.47045)
##  9 Cairns , Australia        17157.  (145.771 -16.92035)
## 10 Canberra , Australia       5953.   (149.13 -35.28094)
## # ... with 39 more rows
YlOrBr <- c("#FFFFD4", "#FED98E", "#FE9929", "#D95F0E", "#993404")

ggmap(map) +
  stat_density_2d(
    data = Rlocations,
    aes(x = lon, y = lat, fill = stat(x)),
    alpha = .1,
    bins = 30,
    geom = "polygon"
  ) +
  scale_fill_gradientn(colors = YlOrBr)

Média de Precipitação por Mês

RD <- dataset[!is.na(dataset$Rainfall), ] %>%
  select(Date, Rainfall)

# substring(RD$Date,6,7) # Month only
RainMonth = aggregate(RD$Rainfall, by=list(Category=substring(RD$Date,6,7)), FUN=mean) # Mean Rainfall by month
ggplot(data=RainMonth, aes(x=Category, y=x, group=1)) +
  geom_line()+
  geom_point()

Criando o modelo

Removendo colunas com muitos dados faltando e depois removendo entradas incompletas.

print(colMeans(is.na(dataset))) # Percentage of missing data in each column
##          Date      Location       MinTemp       MaxTemp      Rainfall 
##    0.00000000    0.00000000    0.01020899    0.00866905    0.02241853 
##   Evaporation      Sunshine   WindGustDir WindGustSpeed    WindDir9am 
##    0.43166506    0.48009762    0.07098859    0.07055548    0.07263853 
##    WindDir3pm  WindSpeed9am  WindSpeed3pm   Humidity9am   Humidity3pm 
##    0.02906641    0.01214767    0.02105046    0.01824557    0.03098446 
##   Pressure9am   Pressure3pm      Cloud9am      Cloud3pm       Temp9am 
##    0.10356799    0.10331363    0.38421559    0.40807095    0.01214767 
##       Temp3pm     RainToday  RainTomorrow 
##    0.02481094    0.02241853    0.02245978
datasetClean <- subset(dataset, select=-c(Evaporation, Sunshine, Cloud9am, Cloud3pm, Date)) # Removing features with high NA precentage
datasetClean <- datasetClean[complete.cases(datasetClean), ] # Removing NAs

One-Hot encoding nas features categóricas.

WD3<-dummy(datasetClean$WindDir3pm)
WD9<-dummy(datasetClean$WindDir9am)
WGD<-dummy(datasetClean$WindGustDir)
Loc<-dummy(datasetClean$Location)

colnames(WD3) <- paste("WindDir3pm", colnames(WD3), sep = "_")
colnames(WD9) <- paste("WindDir9am", colnames(WD9), sep = "_")
colnames(WGD) <- paste("WindGustDir", colnames(WGD), sep = "_")
colnames(Loc) <- paste("Location", colnames(Loc), sep = "_")

datasetD <- cbind(datasetClean, WD3)
datasetD <- cbind(datasetD, WD9)
datasetD <- cbind(datasetD, WGD)
datasetD <- cbind(datasetD, Loc)


datasetD$RainToday <- factor(datasetD$RainToday,
             levels = c('No', 'Yes'),
             labels = c(0, 1))

datasetD$RainTomorrow <- factor(datasetD$RainTomorrow,
                   levels = c('No', 'Yes'),
                   labels = c(0, 1))

df<-tibble(datasetD)

Separando em datasets de treino e teste.

set.seed(1337)
split = sample.split(df$RainTomorrow, SplitRatio = 0.8)
training_set = subset(df, split == TRUE)
test_set = subset(df, split == FALSE)

Treinando o classificador Random Forest.

y <- training_set$RainTomorrow
X <- subset(training_set, select=-c(RainTomorrow))

classifier = randomForest(x = X,
                          y = y,
                          ntree = 100)

classifier
## 
## Call:
##  randomForest(x = X, y = y, ntree = 100) 
##                Type of random forest: classification
##                      Number of trees: 100
## No. of variables tried at each split: 10
## 
##         OOB estimate of  error rate: 14.46%
## Confusion matrix:
##       0     1 class.error
## 0 67205  3120  0.04436545
## 1  9943 10072  0.49677742

Matriz de confusão.

y_pred = predict(classifier, newdata = subset(test_set, select=-c(RainTomorrow)))

# Confusion Matrix
cm = table(test_set$RainTomorrow, y_pred)
print(cm)
##    y_pred
##         0     1
##   0 16830   751
##   1  2413  2591

Calculando a acurácia do modelo.

# Accuracy Score
accuracy = sum(y_pred == test_set$RainTomorrow) / nrow(test_set)
print(accuracy)
## [1] 0.859907

Visualizando as curvas de ROC e Precision e Recall.

precrec_obj <- evalmod(scores =as.integer(y_pred), labels = test_set$RainTomorrow)
autoplot(precrec_obj)